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In recent times increasing numbers of high-rate GPS stations have been installed around the world and 
set-up to provide data in real-time. These networks provide a great opportunity to quickly capture surface 
displacements, which makes them important as potential constituents of earthquake/tsunami monitoring 
and warning systems. The appropriate GPS real-time data analysis with sufficient accuracy for this purpose 
is a main focus of the current GPS research. In this paper we propose an augmented point positioning 
method for GPS based hazard monitoring, which can achieve fast or even instantaneous precise positioning 
without relying on data of a specific reference station. The proposed method overcomes the limitations of 
the currently mostly used GPS processing approaches of relative positioning and global precise point 
positioning. The advantages of the proposed approach are demonstrated by using GPS data, which was 
recorded during the 2011 Tohoku-Oki earthquake in Japan. 

GPS (Global Positioning System) has attracted increasing attention and numerous applications in hazard 
monitoring 14 resulting in a rapid development of high-rate (a 1 Hz) GPS data processing approaches for 
these applications. High-rate GPS observes displacements directly and thus is particularly valuable in 
case of large earthquakes/tsunamis 5 . Consequently in the recent years, dense GPS monitoring networks have been 
built in seismically active regions, e.g., Japan's GEONET (the GPS Earth Observation Network System, http:// 
www.gsi.go.jp/) and UNAVCO's Plate Boundary Observatory (PBO, http://pbo.unavco.org/). These networks are 
complementary to seismic monitoring networks and contribute significantly to earthquake/tsunami early warn- 
ing and hazard risk mitigation 3,6,7 . 

Appropriate and precise GPS real-time data analysis is crucial for the use of the network data for hazard 
monitoring. Currently, the relative baseline/network positioning technique is predominantly used for this pur- 
pose 3,711 . For moderate-to-short baselines, integer ambiguity resolution can be achieved within a few seconds and 
sometimes with only one observational epoch to achieve a high positioning accuracy of a few cm 8 . For the relative 
positioning (RP) technique, GPS data from a network is analyzed simultaneously to estimate station positions. It 
is complicated by the need to assign baselines, overlapping Delaunay triangles 7 , or overlapping sub-networks 12 . 
This is a significant limitation for the challenging simultaneous and precise real-time analysis of GPS data from 
hundreds or thousands of ground stations. Furthermore, intermittent station dropouts complicate the network- 
based relative positioning. Relative positioning also requires a local reference station, which might itself be 
displaced during a large seismic event, resulting in misleading GPS analysis results. The reference station should 
be sufficiently far from the focal region, but must also be part of a sub-network that has relatively short baselines. 
In the case of large earthquakes, such as the Mw 9.0 Tohoku-Oki event at Japan, the reference station may also be 
significantly displaced, even when it is several hundred kilometers away from the event 11 . 

Alternatively, precise point positioning 13 (PPP) can provide "absolute" displacements with respect to a global 
reference frame (defined by the satellite orbits and clocks) using a single GPS receiver. It is more flexible than the 
relative positioning technique and is widely used for hazard monitoring 1417 . However, the PPP method requires a 
long convergence period of about 20 minutes after receiver activation or after serious and/or long signal inter- 
ruption for most of the GPS satellites 18 . The worst case scenario for the GPS component of an earthquake/tsunami 
monitoring system would be a power failure during the disaster, which would reduce the usefulness of the PPP 
based displacement solution because of the time required for re-convergence 19 . To avoid this major disadvantage, 
the PPP regional augmentation 20 has been developed by making use of atmospheric corrections from a regional 
reference network to achieve nearly instantaneous ambiguity resolution. But the regional monitoring stations 
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could also be displaced by the earthquake. Therefore the current PPP 
regional augmentation, in which the reference stations are assumed 
being in static mode and even with known coordinates for generating 
atmospheric corrections 20 or pre-fit undifferenced observation resi- 
duals 25 , could not be used for earthquake monitoring. 

This is our motivation to propose here a novel method for fast or 
even instantaneous positioning, making full use of the currently 
available global PPP service and regional GPS monitoring networks. 
We estimate coordinates of all monitoring stations in kinematic 
mode to avoid the effects of the earthquake induced- displacements 
on atmospheric corrections. The derived atmospheric corrections at 
the stations with fixed ambiguities then can be provided to other 
monitoring stations for instantaneous ambiguity resolution, so that 
precise displacements can always be achieved within a few seconds. 
The series of displacements, derived using the proposed method, will 
be uninterrupted even in case of a break in tracking (loss of signal 
lock, cycle slips, or data gaps) due to a power outage or similar 
disruption. This is a considerable advantage for hazard monitoring 
application. The new method does not depend on a specific reference 
station and therefore the analysis results will not be affected by sim- 
ultaneous shaking of any particular station. It also has better flexibil- 
ity and efficiency compared to complicated network/subnetwork 
analysis. We demonstrate the advantages of the novel augmented 
PPP approach using 1 Hz GEONET data, collected during the 
Tohoku-Oki earthquake (Mw 9.0, 11 March, 2011) in Japan. 



Results 

The 201 1 Mw 9.0 Tohoku-Oki earthquake (11 March 201 1, 05:46:23 
UTC) in Japan is one of the best GPS recorded large earthquakes, as 
Japan has one of the densest GPS networks in the world. The 
Geospatial Information Authority of Japan (GSI) operates more than 
1,200 continuously observing GPS stations (collectively called the 
GPS Earth Observation Network System) all over Japan. The geo- 
graphical distribution of the stations is indicated in Figure 1. The use 
of the GEONET data provides an excellent opportunity to evaluate 
the performance of our novel PPP analysis method. We replayed all 
the 1 Hz GPS data collected by the GEONET stations during the 
2011 Tohoku-Oki earthquake using the augmented PPP method in 
simulated real-time mode. 

First we process 1 Hz GPS ground tracking data of about 80 ~ 90 
globally distributed real-time IGS stations using the GFZ's EPOS-RT 
software in simulated real-time mode for providing GPS orbits, 
clocks and uncalibrated phase delays (UPD) 18,21,23 corrections at a 
5 s sampling interval. Using the orbits, clocks and UPD data, the 
integer ambiguities are fixed in PPP mode for all of the GEONET 
stations and atmospheric corrections are derived on an individual 
station basis. For each GEONET station, three nearby GPS stations 
are selected as augmenting stations. In addition to the orbit, clock 
and UPD data products from the global PPP service, the atmospheric 
corrections of the augmenting stations are interpolated and imposed 
as a constraint on related parameters. Then instantaneous ambiguity 
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Figure 1 | Location of the 2011 Tohoku-Oki earthquake epicenter and the distribution of the high-rate GPS sites. The epicenter is marked by the red 
star. The brown circles represent GPS sites. The black diamond represents the reference site of relative positioning analysis. The purple rectangles 
represent the sites of the time series examples. This figure is drawn using GMT software 32 . 
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(a) Augmented PPP solution at 0176 



(a) Augmented PPP solution at 0175 
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(b) Augmented PPP solution at 0065 
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Figure 2 | The displacement waveforms derived from the augmented PPP 
solution. The north, east and up components are shown by the blue, red 
and black curves, respectively, (a) The displacement waveforms at station 
0176; (b) The displacement waveforms at the reference station 0065. 

fixing is performed independently at each epoch. The displacement 
waveforms, derived from augmented PPP solution, at station 0176 
are shown in the Figure 2a, to illustrate typical behavior. The north, 
east and up components are respectively shown by the blue, red and 
black curves. 

It is found that there are significant GPS data gaps or cycle slips 
during the seismic shaking at some GEONET sites (e.g., 0175, 0588, 
etc). There is for example a data gap of about 2 min at station 0175, 
which starts at epoch 05:47:21 and ends 05:49:27 (GPS time, GPST). 
The displacements from the augmented PPP solution for station 
0175 are shown in Figure 3a. Stations 0172, 0914 and 0918 are 
selected as the augmenting stations for 0175. The estimated iono- 
spheric corrections during seismic shaking at the augmenting sta- 
tions are illustrated in Figure 4a and Figure 5a. The estimated zenith 
wet delays during the 600 s seismic shaking period are respectively 
5.4 ±0.1 cm, 6.5 ±0.1 cm and 6.4 ±0.1 cm at three augmenting 
stations. With the atmospheric corrections, retrieved from the aug- 
menting stations, the atmospheric delays for 0175 are interpolated 
using the linear combination method. The resulting interpolations 
are compared with the estimated values at 0175 in order to assess the 
accuracy of the interpolation. Figure 4b and Figure 5b show the 
ionospheric differences between interpolated and estimated values. 
The differences are found to be smaller than 5 cm. The tropospheric 
interpolation error is about 0.26 cm. We found that the interpolated 
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Figure 3 | Comparisons of the displacement waveforms derived from 
augmented PPP, global PPP and relative positioning solutions. The 

north, east and up components are shown by the blue, red and black 
curves, respectively, (a) The displacements derived from augmented PPP 
solution at station 0175, which has a data gap of about 2 min during the 
seismic shaking; (b) The displacements derived from the global PPP 
solution at station 0175; (c) The displacements derived from relative 
positioning at station 0175, with 0065 as reference station. 



atmospheric corrections are accurate enough for rapid ambiguity 
resolution. 

We also derive displacement waveforms of all GEONET stations 
from the relative positioning (RP) and global PPP solutions, and 
compare them with the augmented PPP solution. The global PPP 
displacements for station 0175 are shown in Figure 3b. In the global 
PPP solution, the displacement series shows a large disturbance after 
the data gap that is caused by the convergence sequence for fixing the 
PPP ambiguities (about 20 min). This unstable behavior is an 
unavoidable problem for a real-time PPP use as the sparse global 
reference network employed cannot provide accurate atmosphere 
delays for fast ambiguity resolution. The relative positioning solution 
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Figure 4 | Ionospheric corrections at augmenting stations and 
ionospheric interpolation errors, (a) The estimated ionospheric 
corrections for GPS satellite PRN 15 at the augmenting stations 0172, 0914 
and 0918 during seismic shaking; (b) Ionospheric interpolation errors of 
PRN 15 for the station 0175. 
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Figure 5 | Ionospheric corrections for the augmenting stations and the 
ionospheric interpolation errors, (a) The estimated ionospheric 
corrections for GPS satellite PRN 28 for the augmenting stations of 0172, 
0914 and 0918 during seismic shaking; (b) Ionospheric interpolation 
errors of PRN 28 for the station 0175. 



for the station 0175 is also shown in Figure 3c. For the relative 
positioning analysis, we adopt the same reference station 0065 as 
Ohta et al (2012) 11 . It can be seen that there are some fluctuations 
in the displacement series derived from the relative positioning solu- 
tion, which are caused by the ground shaking at the reference station 
location. The Figure 2b shows the ground displacements at the 0065 
reference station. Peak surface displacements of up to half a meter 
were recorded at this station during the earthquake even though it is 
about 700 km away from the epicenter. The displacement waveforms 
for the station 0588 derived from augmented PPP (0217, 0590 and 
0965 are selected as augmenting stations), global PPP and RP solu- 
tions are also compared in Figure 6. 

The permanent coseismic displacements of ninety evenly- distrib- 
uted stations derived from post-processed ARIA solution (5 min 
solution), real-time augmented PPP, global PPP, and RP solution 
are shown in Figure 7a, 7b, 7c, and 7d, respectively, by the red arrows. 
The post-processed ARIA solution is provided by the ARIA team at 
JPL (Jet Propulsion Laboratory) and Caltech (California Institute of 
Technology) 29 . It can be found that the permanent coseismic displa- 
cements, derived from the real-time augmented PPP solution, are 
quite consistent with those of post-processed ARIA solution in both 
horizontal and vertical components. The root mean squared errors 
(RMS) of the differences between the two solutions are 1.4, 1.1, and 
1.7 cm in north, east, and vertical components, respectively. 
Figure 7c shows some significant differences between global PPP 



and ARIA displacements at some stations, which are caused by the 
data interruptions at these stations. The corresponding RMS values 
of the differences are 4.3, 22.7, and 9.0 cm in north, east, and vertical 
components. Figure 7d shows, that the RP displacements have obvi- 
ous disagreements with the ARIA results at nearly all stations due to 
problem of the earthquake shaking of the reference station. The RMS 
values of the differences are 10.1, 14.1, and 5.7 cm in north, east, and 
vertical components. Figure 8 shows the displacement differences 
between the ARIA solution and the other three solutions. These 
comparisons show that the augmented PPP method can significantly 
improve the reliability and accuracy of earthquake-induced coseis- 
mic displacements in real-time scenarios. 

We derived four fault slip distributions based on the four different 
GPS analysis techniques introduced above. Identical finite fault para- 
meters are used for the four inversions. Identically as done by Wang 
et al. (20 13) 30 , we employ a slightly curved fault plane, parallel to the 
assumed subduction slab. The dip angle increases linearly from 10° 
on the top (ocean bottom) to 20° at about 80 km depth. To avoid any 
artificial boundary effect, a large enough potential rupture area of 
650 km X 300 km is used. The upper edge of the fault is located 
along the trench east of Japan, on the boundary between the Pacific 
plate and the North American plate. The patch size is about 10 km X 
10 km. The rake angle determining the slip direction at each fault 
patch is allowed to vary between 90° ± 20°. Green's functions are 
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(a) Augmented PPP solution at 0588 
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Figure 6 | Comparisons of displacement waveforms derived from 
augmented PPP, global PPP and RP solutions. The north, east and up 
components are shown by the blue, red and black curves, respectively. 

(a) The displacements derived from augmented PPP solution at station 
0588, which has a data gap of about 2 min during the seismic shaking; 

(b) The displacements derived from global PPP solution at station 0588; 

(c) The displacements derived from relative positioning at station 0588, 
with 0065 as the reference station. 



calculated based on a local CRUST2.0 model by using the software 
codes from Wang et al. (2003) 31 . 

The comparisons of synthetic and observed displacements on 
horizontal and vertical components are shown in Figure 7, and the 



inverted fault slip distributions are shown in Figure 9. Although the 
four results show similar slip distribution, the inversion from real-time 
augmented PPP solution is the most consistent with post-processed 
ARIA solution not only for the slip distribution, but also for the 
displacement fittings. Supposing that the post-processed ARIA result 
is the most reliable and can be taken as a reference for other three 
results, the inversion of global PPP has the worst slip distribution, and 
the inversion of RP solution has the worst displacement fittings. 
Figure 10 shows the fault slip differences between the ARIA solution 
and the other three solutions. Overall, the comparison of the inversions 
shows that the augmented PPP method is beneficial for fault slip 
inversion in real-time scenarios. It provides a more accurate and robust 
estimation of the fault slip distribution and displacement fittings than 
the global PPP solution and RP solution. By contrast, the global PPP 
and RP solutions result in relatively poor slip distributions not only in 
peak slip, but also in the extension of the slip areas (Figure 10). 

Discussion 

We proposed a new GPS analysis method for hazard (e.g. earthquake 
and tsunami) monitoring. The new augmented PPP method can 
overcome the limitations of current relative positioning and global 
PPP approaches for this application. The performance of the new 
approach is evaluated by GPS ground network data, observed during 
the 2011 Tohoku-Oki earthquake in Japan. 

The atmospheric corrections retrieved from the nearby monitor- 
ing stations can be interpolated with accuracy better than 5 cm. This 
means that the interpolated atmospheric corrections are accurate 
enough for rapid ambiguity resolution, which is a prerequisite to 
achieve the most precise displacements. The displacement wave- 
forms, derived using the augmented PPP approach are immune to 
the convergence problem caused by data gaps and cycle slips and the 
problem of the earthquake shaking the reference station compared to 
the waveforms based on RP and global PPP analysis. This makes 
augmented PPP potentially appropriate for the application in opera- 
tional earthquake/tsunami monitoring and warning systems. The 
reliability and accuracy of permanent coseismic displacements are 
also significantly improved. The RMS accuracy of about 1.4, 1.1, and 
1.7 cm are achieved in the north, east, and vertical components, 
respectively. The inversion results indicate that the augmented 
PPP solution is the most consistent with post-processed ARIA solu- 
tion both in the fault slip distribution and displacement fittings. 

Methods 

Successful resolution of integer-cycle carrier-phase ambiguities is a prerequisite to 
achieve the most precise position estimates with GPS by transforming precise but 
ambiguous phase range measurements into precise unambiguous measurements 21,22 . 
For relative positioning, the uncalibrated phase delays (UPD) are removed by the 
application of the double-difference (DD) technique and thus the phase ambiguity 
can be fixed to integers 21 . The atmospheric delays are also mostly eliminated in case of 
moderate-to-short baselines, so that integer-cycle phase ambiguities can be fixed 
within few seconds. Recent studies show that the UPDs can be estimated with high 
accuracy and reliability from a global reference network and transferred to the GPS 
monitoring station to allow resolution of the ambiguities without differencing 18,23 . 
Several international GNSS service (IGS) analysis centers provide GPS orbit, clock, 
and UPD data products to allow real-time PPP use enabling ambiguity resolution 
anywhere in the world 24 " 26 . However, PPP still needs a comparatively long (re)con- 
vergence time of approximate 20 minutes to achieve reliable integer ambiguity 
resolution because precise atmospheric delay models cannot be derived from such a 
sparse global reference network 18 . 

An increasing number of regional GPS monitoring networks are installed around 
the world for precise navigation and geophysical applications, especially in seismic- 
ally active regions (e.g. fapan, Western North America, Greece, and Chile). One 
possible solution for achieving fast ambiguity resolution in PPP is to retrieve the 
atmospheric delays as corrections from data of these dense regional networks. By 
applying the UPD corrections, the integer un-differenced ambiguities on the LI and 
L2 frequencies can be fixed in PPP mode at all regional monitoring stations. The 
atmospheric corrections of the ionospheric slant and tropospheric zenith wet delay 
then can be derived from the PPP fixed solution as 
I' r ,-m s r -Z r + u s r -Ar= -f r 



- f + t r + Xft> rj - bf) + XjN s rj + $ j 



(1) 
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Figure 8 | The residual displacements from the ARIA solution, (a) Residual differences between augmented PPP and ARIA vectors; (b) Residual 
differences between global PPP and ARIA vectors; (c) Residual differences between RP and ARIA vectors. This figure is drawn using GMT software 32 . 



- I'j -m s r -Z r + u* r -Ar= -p s rj -f + t r + c(d rJ + df) + f? rj 



(2) 



Where, l s r j, p s r , denote phase and code observables from satellite s to receiver r at 
frequency j; u\ is the unit direction vector from site to satellite; Ar denotes the 
increments of the receiver positions; Z r denotes the tropospheric zenith wet delay; m s r 
is the wet mapping function; f and t r are the clock errors; Xi is the wavelength; b r i is the 
receiver-dependent uncalibrated phase delay; bf is the satellite- dependent UPD; d r i 
and df are the code biases; I s r j is the ionospheric delay; N*j is the integer phase 



ambiguity; e s r - and are the measurement noise terms of the pseudo-range and 
carrier phase. 

This procedure is very flexible and computational efficient to be applied even for 
monitoring networks with a large number of stations as the atmospheric corrections 
are derived for each station individually. Because regional monitoring stations 
themselves could be displaced by the earthquake, the coordinates are estimated in 
kinematic mode to avoid the effects of earthquake induced -displacements on the 
atmospheric corrections that are generated. The constraints imposed on the 
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kinematic coordinates of adjacent epochs are fine-tuned by using an adaptive filter 27 
in real-time to strengthen the solution. Usually atmospheric delay is rather stable over 
short periods and can be represented by a constant or a linear function. Therefore, 
even in periods of strong shaking, station position and atmosphere are distinguish- 
able in parameter estimation because of the significant difference in their temporal 
characters. 

A polynomial model can be used to represent the derived atmospheric corrections 
on small regional scales. Here three or more nearby monitoring stations are selected 
as augmenting stations for each monitoring station, and the atmospheric corrections 
of the selected augmenting stations are interpolated by a Linear Combination Method 
(LCM) 28 as 

V) «, = 1, ^ <*i(Xm -Xf) = 0, V) y-i 2 =Min (3) 
;=i i=i i=i 



v m =Y; a iVi (4) 

Where, n is the number of selected augmenting stations; m and /' are indices for the 
monitoring and the selected augmenting stations, respectively; a, denotes the inter- 
polation coefficient; X m and X, are the station coordinates in the local horizontal 
plane system; AX im and AY im are the plane coordinate differences between the 
monitoring and augmenting station; v, is the ionospheric or tropospheric delay; v m is 
the interpolated ionospheric or tropospheric delay at the monitoring station. 

For regional reference networks with moderate-to-short baselines (few tens of 
kilometers inter-station distance) cm-level accuracy can be achieved for the inter- 
polated atmospheric delay corrections. The precise interpolated atmospheric cor- 
rections are imposed as a strong constraint on the related parameters of the 
monitoring station, while the coordinates are estimated in kinematic mode. 
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Figure 10 | The residual slip distributions from the ARIA solution, (a) Residual differences between augmented PPP and ARIA inversion; (b) Residual 
differences between global PPP and ARIA inversion; (c) Residual differences between RP and ARIA inversion. The star denotes the epicenter. This figure is 
drawn using GMT software 32 . 



Assuming that r x to r n are selected as augmenting stations for interpolating correc- 
tions for the monitoring station r m . The ionospheric slant delay parameter for an 
individual satellite s> is constrained to the interpolated correction as 

?„ m~N(0,<) (5) 
And the constraint for the zenith wet delay parameter is 

Zr m -Z ri ,ri-m = w r, wt~N(0,<j 2 Wt ) (6) 

Where Jp denotes the slant ionospheric delay from station r m to satellite 5,; I* ri _ r is 
the interpolated ionospheric correction; Z r denotes the zenith wet delay for station 
r m , and Z r] fa ... r is the interpolated correction. Wj and Wt are the biases between the 
true and the interpolated atmospheric corrections. The statistical processes of Wj and 
Wt are zero mean white processes with variance of o~ z y} and o 2 Wr for the ionospheric 
and tropospheric delays, respectively. 

By adding this precise atmospheric delay model to the orbit, clock and UPD 
products used in global PPP ambiguity resolution, instantaneous ambiguity resolu- 
tion is achievable at the monitoring station, so that the augmented PPP can have 
ambiguity resolution performance equivalent to relative positioning. It should be 
mentioned that the selection of augmenting stations is critical, as atmospheric cor- 
rections can only be derived from augmenting stations at which the ambiguity 
resolution is successfully achieved. 
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